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1.0  INTRODUCTION 


1.1  Objective 

Because  Computational  Fluid  Dynamics  (CFD)  can  treat  the  nonlinear  phenomenon 
in  high-Reynolds-number  flow  simulation  with  the  acid  of  high  speed  digital  computer, 
CFD  are  becoming  in  the  matured  stage  for  both  physical  research  problems  and  practical 
engineering  problems,  it  is  becoming  an  integral  part  of  the  study  of  fluid  mechanics  and 
the  aerodynamics  design  process  [39].  Research  topic  on  CFD  are  generally  directed 
toward  solving  problems  for  understanding  flow  physics,  treating  realistic  configurations, 
maintaining  good  computational  efficiency  and  stability,  postprocessing  of  the  computer 
results,  and  maintaining  user  friendliness.  Many  researchers  have  been  searching  for 
methods  to  simulate  flows  around  very  complex  geometries,  such  as  the  integrated  space 
shuttle  configuration  [Ref.  21],  and  full  aircraft  configurations  [Ref.  25],  which  depends 
on  the  use  of  different  technologies  developed  for  CFD  and  other  relative  fields,  such  as 
grid  generation,  flow  solver,  and  flow  visualization. 

Currently,  with  limited  computation  resources,  the  methods  of  treating  realistic 
configurations  and  increasing  accuracy  of  the  numerical  simulation  are  one  of  the  most 
important  tasks  in  the  computational  fluid  dynamics  research  and  practical  engineering 
application.  For  a  single  body  shape,  as  shown  in  Figure  1  for  space  shuttle  orbiter,  the 
use  of  a  single  body-conforming  curvilinear  mesh  leads  to  the  most  efficient  solution 
procedure.  Flowever,  for  more  complex  shapes,  especially  for  complex  three-dimensional 
configuration,  as  shown  in  Figure  2,  for  a  simplified  space  shuttle  configuration,  it  is  quite 
difficult  to  generate  a  single  body-conforming  well-clustered  curvilinear  grid  that  is  not 
overly  skewed  and  has  smooth  variation.  On  the  other  hand,  the  accuracy  enhancement  is 
another  obstacle  in  the  computational  fluid  dynamics.  No  matter  what  discretization 
methods  would  be  used,  finite  difference  or  finite  volume,  the  solution  of  the  partial 


1 


Figure  2.  Complex  space  shuttle 
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differential  equations  of  fluid  motion  requires  that  the  computational  domain  and  dependent 
variables  be  represented  on  a  network  of  discrete  points.  The  distribution  of  these  points  is 
influenced  by  the  choice  of  the  coordinate  system,  the  order  of  the  numerical 
approximation,  and  the  location  of  strong  geometric  and  flow  field  gradients.  As  long  as 
the  discretization  method  is  used,  the  solutions  of  the  differential  equations  become  more 
accurate  as  the  number  of  the  grid  points  are  increased.  But  because  of  the  limitation  of 
computer  resources,  it  is  necessary  to  improve  the  grid  resolution  without  increasing  the 
number  of  grid  points. 

Many  different  methods  have  been  developed  in  grid  generation  and  discretization 
to  improve  the  grid  resolution.  Grid  adaptation  is  one  of  the  solutions.  Local  grid 
resolution  can  be  improved  by  adopting  solution-adaptive  grid  technique  without  increasing 
the  total  number  of  grid  points.  Another  method  is  local  grid  enrichment  In  addition  to  the 
original  grid,  locally  defined  finer  grid  is  created  and  overlapped  or  overlaid  on  the  region 
of  interest  of  the  original  grid. 

For  geometry  imposed  complex  flow  fields,  composite  strucmred  grid  schemes  and 
unstructured  grid  schemes  are  two  major  approaches  for  nonlinear  simulation.  For  the 
most  part,  the  unstructured  grid  approaches  have  come  from  the  finite  element  community, 
although  finite  difference  and  finite  volume  methods  can  be  written  for  an  unstructured 
grid.  Both  approaches  require  a  more  complex  data-handling  program  than  that  required  by 
a  simple  single-structured  grid,  and  both  approaches  have  their  strengths  and  weaknesses. 
It  seems  that  unstructured  grid  methods  are  generally  considered  to  be  more  versatile  and 
easier  to  adapt  to  complex  geometry,  but  it  has  difficulty  in  treating  viscous  flow  problems, 
while  composite  structured  grid  methods  are  generally  considered  to  use  more  efficient 
numerical  algorithms,  require  less  computer  memory,  and  can  handle  viscous  flow 
problems  without  difficulty.  The  best  approach  is  still  not  resolved,  but  hybrid  schemes, 
which  use  is  structured  grid  in  most  areas  and  only  unstructured  grid  in  limited  regions, 
have  already  appeared  [Refs.  40  through  42],  it  incorporates  the  best  features  of  both 
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schemes,  and  will  certainly  grow  in  popularity  and  ultimately  prove  to  be  the  best  for 
treating  complex  configurations. 

Composite  grids  use  more  than  one  grid  to  mesh  an  overall  configuration  with  each 
individual  grid  of  the  system  patched  or  overset  together.  Patched  grids  are  individual 
meshes  that  are  joined  together  at  some  common  interface  plane,  but  it  is  still  relatively 
difficult  in  grid  generation  because  various  interfaces  have  to  be  defined  and  grids  have  to 
be  generated  with  both  inner  and  outer  defined  boundary  surfaces.  Overset  grids  do  not 
require  a  common  interfaces,  but  rather,  a  simply  superimposed  or  partially  superimposed 
mesh  is  needed  to  cover  the  region  of  interest,  and  to  provide  the  means  of  matching  the 
solutions  across  boundary  interface.  Without  any  special  way  of  grid  connection  or 
common  interface,  the  grid  generation  is  greatly  simplified  since  the  hyperbolic  grid 
generation  technique  may  be  used. 

The  composite  overset  or  Chimera  method  can  provide  the  flexibility  to  employ 
boundary  conforming  grids  on  component  parts  of  the  geometry,  define  the  mesh 
selectively  in  the  regions  of  interest,  and  permit  different  flow  models  to  be  solved  on  each 
grid.  In  this  discretization  process,  each  grid  in  the  system  is  ordered  and  is  thus  suitable 
for  efficient  finite  difference  solutions  on  vector  computer  and  any  single-grid  code.  The 
primary  objective  of  present  study  is  to  develop  a  composite  overset  or  Chimera  code 
which  can  handle  general  3-D  configuration  problems. 

For  a  complex  geometric  configuration  problem,  to  generate  composite  overset 
grids,  and  construct  required  data  structure  for  Chimera  scheme  is  not  an  easy  job,  and  it  is 
more  costly  to  run  a  three-dimensional  Chimera  scheme  for  Navier-Stokes  equations. 
Therefore,  before  running  the  Navier-Stokes  solver,  it  is  very  important  to  check  Chimera 
data  structure  to  make  sure  that  all  the  Chimera  overset  grids  connections  are  working 
correctly.  Therefore,  besides  developing  a  Chimera  connecting  code,  the  other  purpose  of 
present  work  is  to  develop  a  simple  numerical  procedure  in  three  dimension,  which  can 
easily  fit  in  Chimera  scheme  and  give  an  initial  inviscid  flow  solution.  These  inviscid 
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results  from  the  simple  flow  solver  can  be  used  to  check  the  Chimera  scheme,  and  to  be 
used  as  the  initial  results  of  Navier- Stokes  solver. 

1.2  Literature  Review 

With  the  aid  of  supercomputer,  CFD  researchers  start  to  calculate  complex 
configurations,  e.g.,  Meakin,  et  al.  [Refs.  17,  18,  and  19],  used  overset  grid  method  to 
simulate  transient  flow  field  by  solving  Navier-Stokes  equations  for  integrated  space  shuttle 
vehicle  configuration  during  ascent  and  SRB  separation.  Buning,  et  al.  [Refs.  21,  22  and 
23],  used  an  overset  grid  method  to  solve  the  thin-layer  Navier-Stokes  equations  for  the 
integrated  space  vehicle  launching  configuration.  Flores,  et  al.  [Ref.  56],  solved  the 
Navier-Stokes  equations  for  a  transonic  flow  by  using  a  patched  grid  method.  Obayashi 
[Ref.  57]  solved  Navier-Stokes  equations  for  a  wing-fuselage-tail  geometric  configuration 
by  using  a  single  grid  system,  and  Jameson,  et  al.  [Ref.  58],  solved  Euler  equations  for 
Boeing  747-200  using  an  unstructured  grid. 

Of  these  calculations  of  complex  geometries,  the  grid  system  can  be  roughly 
divided  into  the  following  types: 

•  single  grid 

•  composite  grids,  which  can  be  divided  into: 

-  patched  grid 

-  overset  grid 

-  adaptive  grid 

For  a  single  grid,  it  is  simple  to  solve  control  equations  but  it  is  difficult  for  grid 
generation  for  a  complex  geometric  configurations,  and  it  may  contain  overly  skewed 
meshes,  which  could  give  an  inaccurate  solution.  On  the  other  hand,  using  a  single  grid 
could  not  refine  the  grid  selectively,  therefore,  it  may  require  more  grid  points  than  a 
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comparable  multiple  grid  approach.  Comparing  with  a  single  grid,  a  multiple  grid  method 
is  flexible,  but  it  also  demands  a  significant  amount  of  computer  resources  and  different 
techniques,  which  may  involve  computational  domain-decomposition,  connection  or 
communication  between  each  grid  and  composite  grids  data  management 

Domain-decomposition  techniques  subdivide  the  computational  domain  into 
simpler  subdomains  which  admit  a  more  easily  constructed  mesh.  Several  strategies  have 
been  explored  to  subdivide  the  domain  and  establish  communications  among  the 
subdomains.  One  group  of  approaches,  the  grid-patching  or  zonal  methods,  uses  common 
or  shared  boundaries  and  another  uses  embedded  or  overset  grids  to  subdivide  the  domain. 
The  work  of  Rubbert  and  Lee  [Ref.  56]  is  typical  of  the  methods  which  construct  a  global 
mesh  from  subdomains  which  share  common  subdomains.  They  generate  a  global  mesh 
by  solving  grid-generation  equations  on  all  subdomains  simultaneously  and  by  requiring 
that  the  grid  lines  be  continuous  across  subdomain  boundaries.  A  difficulty  with  this 
approach  is  that  irregularities  which  occur  in  comers  and  along  boundaries  impose 
constraints  on  the  algorithm  used  to  solve  the  flow  equations.  Lasinski  [Ref.  43]  takes  can 
alternate  approach  and  solve  for  the  flow  field  on  each  subdomain  separately  with 
communication  among  the  grids  established  by  the  transfer  of  boundary  data.  In  their 
approach,  the  patches  overlap  one  point  with  common  points  on  the  boundary  to  obviate 
interpolation  for  boundary  data.  Hessenius  and  Pulliam  [Ref.  44]  have  modified  the 
approach.  Rai  [Ref.  45]  further  generalized  the  method  to  admit  independent  grids  in  each 
subdomain.  Communication  across  grid  boundaries  is  accomplished  by  means  of  special 
difference  formula  at  the  boundaries  which  maintain  conservation  properties  across  the 
subdomains.  Similar  methods  have  been  developed  by  Miki  and  Takagi  [Ref.  46}.  Holst, 
et  al.  [Ref.  47]  have  applied  the  technique  to  large  3-D  grids. 
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The  grid-embedding  or  oversetting  techniques  do  not  require  common  boundaries 
between  subdomains,  but  rather,  a  common  or  overlap  region  is  required  to  provide  the 
means  of  matching  the  solutions  across  boundary  interfaces.  The  usual  procedure  uses 
interpolations  of  embedded  boundaries  to  provide  the  necessary  communication  among  the 
grids.  There  are  several  implementations  of  the  method.  Steger,  et  al.  [Ref.  33]  and 
Benek,  et  al.  [Ref.  4],  developed  a  "Chimera"  scheme  in  three  dimensions  for  the  solution 
of  both  a  linearized  flow  model  or  Euler  equations.  Dougherty  [Ref.  38]  has  extended  the 
grid-embedding  technique  to  time-dependent  motions.  Atta  {ref.  48]  applied  this  method  to 
solve  the  full-potential  equation  in  two  and  three  dimensions  using  a  separate  implicit 
solution  algorithm  for  each  grid.  Fuchs  [Ref.  49]  employed  this  approach  to  internal 
flows,  and  Rai  (Ref.  50}  used  a  combination  of  patched  and  overset  grids  to  solve  rotor- 
stator  interaction  problems. 

The  second  technique,  grid-adapting  methods,  forces  the  mesh  to  evolve  with  the 
solution  of  the  flow  equations,  the  advantage  of  this  method  is  that  the  initial  mesh  does 
not  need  to  anticipate  accurately  all  regions  of  large  flow  gradients,  thus  it  can  make  the 
most  efficient  use  of  available  mesh  points,  and  reduce  the  grid-generation  effort  by 
automatically  clustering  grid  points  to  regions  of  high  gradients.  Gnoffo  [Ref.  51]  models 
the  mesh  as  a  network  of  springs  whose  constants  are  determined  from  the  flow  gradients. 
Nakahaski,  et  al.  [Ref.  52],  extend  these  methods  to  both  steady  and  unsteady  flow 
problems.  Ghia,  et  al.  [Ref.  53],  used  both  the  grid-evolution  equation  and  the  flow 
equation  by  requiring  that  the  coefficient  of  the  convective  term  in  the  flow  model  be 
minimized.  Brackbill  and  Saltzman  [Ref.  54]  used  variational  methods  to  produce  grid- 
evolution  equations.  Berger  [Ref.  56]  developed  a  dynamic  grid  refinement  technique  by 
successively  embedding  finer  grids  to  resolve  flow  gradients. 
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2.0  CHIMERA  APPROACH  AND  GRID  TOPOLOGY 


2.1  Introduction 

For  complex  geometries,  grid  generating  itself  is  not  an  easy  task.  Though  it  is 
possible  to  generate  a  single  grid  for  a  complex  geometry,  the  resultant  grid  is  most  often 
overly  skewed  in  certain  directions  or  regions,  or  doesn't  have  the  needed,  clustering  to 
resolve  the  flow  field  in  regions  of  rapid  change.  One  good  method  to  overcome  this 
difficulty  is  to  divide  the  complex  shape  into  several  simple  ones  and  generate  the  grid 
about  these  simple  shapes,  then  overset  them  one  on  the  top  of  the  others.  With  using 
inter-grid  communication  in  the  flow  solver,  the  overset  grid  can  be  used  to  solve  complex 
geometric  problems. 

The  overset  grid  approach  used  in  this  research  was  first  devised  by  Steger,  and  his 
colleague  [Refs.  2  through  5,  10,  16,  21,  22,  and  33,  etc.],  and  given  the  name,  the 
Chimera  approach,  after  the  Greek  legendary  creature  that  was  compounded  of 
incompatible  parts  which  signifies  that  the  Chimera  approach  can  take  incompatible  grids 
(i.e.,  no  common  boundaries  between  different  grids)  and  'glue'  them  together  to  be  solved 
by  the  flow  solver.  Because  the  Chimera  approach  does  not  require  common  boundaries 
between  component  meshes,  complex  grid  generating  can  be  replaced  by  a  combination  of 
a  series  of  simple  ones.  Therefore,  it  provides  the  flexibility  to  employ  boundary- 
conforming  grids  on  component  parts  of  the  geometry,  to  change  grids  of  some  component 
without  changing  other  component  grids.  Furthermore,  it  is  also  easy  to  refine  the  mesh 
selectively  in  the  regions  of  interest  without  regenerating  the  whole  grid  for  the  complete 
configuration.  Another  feature  of  the  Chimera  approach  is  that  the  flow  solver  simulation 
is  done  in  sequence  for  component  grids,  which  offers  a  saving  in  memory  usage  for 
solving  flows  around  a  complex  geometry  since  it  only  requires  memory  enough  to  handle 
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the  largest  component  grid,  and  it  is  also  possible  to  use  different  computers  to  update  the 
results  on  each  grid  if  a  multitasking  computer  processor  is  available. 

In  this  discretization  process,  each  individual  grid  in  the  system  is  well  ordered  and 
thus  suitable  for  efficient  finite  difference  solution  using  highly  vetorized  computers  and 
any  available  single-grid  code.  It  is  also  possible  to  use  different  models  or  schemes  with 
different  time  steps,  etc.,  for  different  components,  which  opens  many  possibilities  for 
enhancing  the  rate  of  convergence. 

The  Chimera  method  has  been  demonstrated  for  a  wide  variety  of  aerodynamic 
applications  including,  among  many  others,  simulations  about  the  integrated  space  shuttle 
vehicle  for  various  ascent  conditions,  [Refs.  21  through  24],  transonic  computations  about 
the  entire  F-15  aircraft  with  pylongs  and  stores  {Ref.  12],  simulation  of  the  viscous  flow 
about  the  F/A-18  at  high  angles  of  attack  [Ref.  27],  and  numerical  studies  of  a  delta  wing 
platform  with  multiple  jets  in  ground  effects  [Ref.  7].  In  addition,  overset  grid  techniques 
represent  a  powerful  tool  for  analyzing  problems  involving  relative  motion  between  vehicle 
components.  Such  applications  have  been  carried  out  time-accurately  in  three  dimensions 
for  the  separation  sequence  of  the  space  shuttle  solid  rocket  boosters  [Ref.  16.17  and  19], 
and  aircraft  store  separation  sequences  [Ref.  10, 1 1  and  19].  The  Chimera  method  has  also 
been  successfully  applied  to  many  nonaerodynamic  problems,  such  as  applications  in 
biomedical  fluid  mechanics  [Ref.  15],  environmental  flow  simulations  [Ref.  18],  and  the 
high  speed  train  project  in  Japan  [Ref.  13]. 

Although  there  are  many  advantages  to  the  overset  mesh  method,  there  are  also 
disadvantages  which  may  influence  the  method  widely  use.  Because  interpolation  points 
and  blanked  points  must  be  located  and  labeled  for  special  treatment,  a  massive  data  base  is 
needed.  These  schemes  are  more  complex  relative  to  using  a  single  grid  scheme  and  data 
interpolation  also  can  cause  inaccuracies  such  as  local  loss  of  numerical  conservation  of 
fluxes.  For  unsteady  problems  with  relative  motion  between  grids,  the  interpolation  data 
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and  index  need  to  be  calculated  at  each  time  step,  therefore,  it  would  be  more  expensive 
than  a  single  grid  scheme. 

Although  there  are  many  applications  of  the  Chimera  method,  the  successful  code  is 
still  very  limited,  especially  for  viscous  flow  problems.  A  very  good  code  PEGSUS 
(versions  2.0,  3.0  and  4.0),  which  was  developed  by  CALSPAN  of  AEDC,  has  been 
widely  used  by  Chimera  group  researchers,  and  has  many  successful  applications  [Refs. 
21  through  23],  but  it  still  has  some  weaknesses  needed  to  overcome,  for  example,  as 
shown  in  Figure  3,  the  space  above  a  concave  body  surface  is  covered  by  two  overset 
grids.  The  values  on  the  boundary  of  the  solid  line  grid  are  interpolated  from  the  solution 
of  the  dashed  line  grid,  and  the  values  on  the  boundary  of  the  dashed  grid  are  interpolated 
from  the  solution  of  the  solid  line  grid.  For  PEGSUS  code,  the  points  below  the  surface 
grid  are  interpolated  from  the  solution  of  solid  line  grid.  For  PEGSUS  code,  the  points 
below  the  surface  grid  are  considered  inside  the  surface,  so  these  points  are  blanked  out. 
For  viscous  flow  problems,  the  grids  are  very  fine  near  the  surface,  thus,  up  to  10  points 
below  the  surface  grid  could  be  blanked  out.  No  information  can  be  obtained  for  these 
points.  The  other  limitation  is  for  high  Reynolds  number  viscous  flow  above  body 
surface,  as  shown  in  Figure  4,  in  the  region  very  close  to  the  surface,  especially  in  the 
sublayer  region  of  the  boundary  layer,  regular  tri-linear  interpolation  method  cannot  give 
satisfied  results,  because  the  value  of  boundary  point  of  the  dashed  line  grid  is  interpolated 
from  the  solutions  of  the  solid  line  grids,  which  are  not  at  the  same  distance  from  the 
surface.  This  is  because  for  high  Reynolds  number  problems,  the  velocity  profiles  have 
very  high  gradient  in  the  surface  normal  direction;  a  small  difference  in  the  surface  normal 
direction  can  bring  out  a  big  difference  in  velocity  value,  so  the  normal  interpolation 
method  needs  to  be  modified.  Therefore,  it  is  necessary  to  develop  new,  more  efficient 
methods  and  logic,  which  can  handle  more  general  and  complex  problems,  such  as  viscous 
flow  problems. 
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Figure  3.  Concave  surface  grids 


Figure  4.  High  Reynolds  number  boundary  layer 
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2.2  Chimera  Method 

The  grid  embedding  has  two  major  parts:  decomposition  of  the  domain  into 
overlapping  subdomains  and  communication  among  the  subdomains.  Each  subdomain  is 
chosen  to  lessen  the  effort  of  generating  an  acceptable  mesh,  each  of  these  meshes 
represents  a  separate  and  independent  grid  generation.  To  increase  the  flexibility  in  the 
domain  decomposition,  it  is  necessary  to  remove  the  regions  of  a  mesh  containing  an 
embedded  grid  from  that  mesh,  that  is,  an  embedded  mesh  introduces  a  'hole'  into  the 
mesh  in  which  it  is  embedded.  Because  appropriate  boundary  values  are  interpolated  from 
the  mesh  in  which  the  boundary  is  embedded,  and  these  parts  do  not  enter  into  the  solution 
process,  inter-grid  communication  is  simplified  through  the  grid  boundaries.  Therefore, 
the  Chimera  procedure  naturally  separates  into  two  parts:  (1)  generation  of  the  composite 
mesh  and  associated  interpolation  data  and  (2)  solution  of  the  flow  equations  on  the 
composite  and  interpolation  data  mesh.  Each  part  is  embodied  in  a  separate  computer  code. 
The  first  code  takes  the  independently  generated  composite  mesh  as  input  and  automatically 
constructs  the  composite  mesh  as  an  output,  which  requires  establishment  of  the  proper 
lines  of  communication  among  the  grids  through  tri-linear  interpolations  within  the 
overlapping  regions  to  determine  outer  boundary  point  properties  for  each  grid.  In  addition 
to  the  outer  boundary  point  communication,  holes  arise  whenever  the  solid  surfaces 
contained  by  one  grid  are  also  contained  within  other  grids,  and  the  hole  requires  boundary 
condition;  therefore,  it  is  necessary  to  identify  points  and  locate  points  from  which 
boundary  values  can  be  interpolated,  then  evaluate  interpolation  parameters.  Because  the 
hole  point  and  hole  boundaries  point  are  not  at  the  computation  domain,  they  are  imposed 
by  blanking  out  the  control  equations  within  the  solution  matrix  and  directly  substituting  the 
boundary  condition.  A  detailed  discussion  of  the  methodologies  used  to  tag  hole  boundary 
points  and  outer  boundary  points  is  given  in  following  sections  2.2  and  2.3. 
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Once  the  subgrids  were  constructed,  and  the  interpolation  coefficients  were 
computed  and  stored  by  the  assembly  code  described  above,  the  second  code  can  use  the 
stored  interpolation  coefficients  to  provide  boundary  condition  for  each  grid,  and  solve  the 
appropriate  flow  equation  on  each  grid  with  single  grid  solver.  A  detailed  discussion  is 
given  in  Section  2.4. 

2.2.1  Hole  Creation 

Because  each  component  mesh  is  generated  independently,  complications 
frequently  arise  when  the  grids  are  embedded.  For  example,  points  of  an  enclosing  mesh 
G,  j ,  may  be  found  to  lie  within  a  solid  boundary  contained  within  an  embedded  grid 
G|+i  i  ,  such  points  lie  within  the  solid  boundary  and  outside  the  computational  domain, 
therefore,  must  be  excluded  from  the  solution  process.  Thus,  the  embedded  mesh  G|^.i ; 
introduces  a  "hole"  in  G,  j.  The  only  computational  requirement  is  that  there  remain  a 
sufficient  overlap  to  support  an  interpolation  for  the  outer  boundary  of  G,+i j  from  points  in 

Gu- 

There  is  a  two-step  procedure  in  finding  the  hole  points  as  well  as  the  hole 
boundary.  First,  (1)  find  a  geometric  center  of  embedded  body  surface  as  a  temporary 
origin,  say  Pq  (2)  find  the  Rmax  and  Rmin  as  the  maximum  and  minimum  distances  from 
Po  to  points  on  embedded  body  surface.  All  the  points  that  fall  outside  of  the  larger  sphere 
with  radial  of  Rmax  are  considered  to  be  field  points.  All  the  points  that  fall  inside  of  a 
smaller  sphere  with  radial  of  Rmin  are  considered  to  be  "hole"  points,  those  points  that  fall 
between  the  two  spheres  need  a  more  accurate  method  to  tell  whether  or  not  a  point  is  a 
hole  point.  To  avoid  unnecessary  confusion  in  a  three-dimensional  case,  a  two- 
dimensional  case  is  presented  to  illustrate  the  idea  of  the  hole  searching  procedure.  In 
Figure  5,  construct  outward  normal  to  the  initial  hole  boundary  at  each  point  defining  the 
surface,  C.  Test  the  magnitude  of  r,  the  position  vector  relative  to  P©.  For  every  point  P 
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Figure  6.  Construction  of  fringe  points 
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of  G,  i ;  if  the  |r|  <  .  P  niay  be  inside  C.  So,  form  N  •  Rp  where  N  is  the  unit 

normal  at  the  surface  point,  Pc,  the  closet  point  to  P,  and  Rp  is  the  position  vector  of  P 
relative  to  Pc;  if  N-Rp  <  0,  P  lies  within  the  hole  and  the  flag  variable,  IB  LANK,  is  set  to 

zero,  otherwise,  P  is  outside  the  hole  and  IB  LANK  is  set  to  1.  The  IBLANK  variable  is 
used  by  the  flow  solver  to  determine  whether  a  point  should  enter  into  the  solution  process 
or  not,  as  illustrated  by  the  equation  below: 

AAQ  =  IBLANK  *  RHS 

where  A  is  the  coefficient  matrix,  AQ  is  the  change  in  flow  solution,  and  RHS  is  the 
source  term.  Thus  for  the  hole  points,  the  above  equation  reduces  to 

AQ=0 

and  the  values  of  the  variable  at  the  hole  points  are  not  changed  in  the  solution  process. 

The  points  of  G,  i  within  a  hole  are  excluded  from  the  solution  and  are  not  usable  as 
boundary  points.  Therefore,  additional  points  of  G,  j  are  identified  as  hole  boundary  or 
fringe  points.  After  the  hole  points  are  found,  the  hole  boundary  points  or  fringe  points 
can  easily  be  located  by  searching  the  IBLANK  values  of  neighboring  points.  Figure  6 
illustrates  the  construction  of  fringe  points.  The  procedure  is  to  examine  the  search  stencil 
for  each  point,  P,  in  G, ^  at  which  IBLANK=1;  if  one  of  the  nearest  neighbors  is  a  hole 

points,  P  is  a  fringe  point.  The  indices  of  the  fringe  points  are  added  to  a  list  of  boundary 
points  which  require  interpolation  data.  The  values  of  the  unknowns  at  these  boundary 
points  are  interpolated  from  the  embedded  mesh,  Gj  j+j  . 

2.2,2  Searching  for  Interpolation  Elements 

Communication  among  the  grid  is  achieved  by  interpolation  of  boundary  values 
from  the  grids  in  which  the  boundaries  are  contained.  Care  must  be  taken  in  the  application 
of  the  method  to  ensure  that  the  interpolated  point  lies  within  this  cell. 
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Figure  7.  A  sphere  surface  covered  by  a  single  piece  of  grid. 
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Figure  8.  The  whole  sphere  surface  covered  by  two  pieces  of  grids. 
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Figure  9.  The  whole  sphere  surface  covered  by  four  pieces  of  grids. 
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For  the  overset  grid  schemes  interpolation  is  generally  needed  to  update  boundary 
points,  and  interpolation  requires  routines  to  seek  out  nearest  points  (or  cells)  to  interpolate. 
Finding  these  points  can  be  costly  and  prone  to  error.  Moreover,  as  Figures  7,  8  and  9 
indicate,  even  for  the  nearest  points  in  space,  it  is  not  necessarily  that  the  best  interpolation 
results  are  obtained. 

There  are  two  simple  methods  which  can  be  used  to  find  the  cell  which  contains 
given  boundary  grid.  One  of  the  methods  is  developed  in  the  physical  coordinate  x,y,z 
space,  which  is  called  "stencil"  method.  The  other  method  is  developed  in  the 
computational  coordinate  T|,  ^  space,  which  is  called  Jacobin  matrix  method,  because  it 
needs  to  use  Jacobin  matrix  to  calculate  the  interpolation  coefficients.  For  our  case,  the 
Jacobin  matrix  method  is  more  efficient  than  "stencil"  method,  because  the  Jacobin  matrix 
was  calculated  in  other  subroutines  in  the  same  program. 

For  more  than  two  mesh  problems,  it  is  possible  that  the  boundary  of  each  mesh  is 
overlapped  by  other  two  or  more  meshes,  which  means  for  the  same  point  of  a  given 
boundary,  its  values  can  be  interpolated  from  different  meshes,  therefore,  it  is  necessary  to 
introduce  new  logic  and  a  series  of  parameters  to  judge  for  which  grids  the  interpolated 
boundary  data  are  best;  therefore,  the  data  structure  becomes  very  complicated.  On  the 
other  hand,  the  simple  logic  can  fail  because  of  the  need  to  introduce  a  starting  point,  which 
is  not  far  from  the  destination  point  and  in  the  same  grid.  For  the  general  multiple  overset 
grids,  these  conditions  are  not  always  satisfied. 

2.2.3  Interpolation  Points 

Because  the  separate  grids  are  to  be  treated  as  independent  entities,  boundary 
conditions  must  be  supplied  to  each.  The  boundary  conditions  of  the  differential  equations 
which  model  the  flow  provide  data  only  at  the  boundaries  of  the  computational  domain. 
Thus,  other  data  must  be  obtained  for  the  subdomain  boundaries  which  are  not  coincident 
with  those  of  the  computational  domain.  Because  the  subdomain  boundaries  typically  lie  in 
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the  interior  of  the  computational  domain  where  the  differential  equations  are  valid,  it  seems 
appropriate  that  the  solution  of  these  equations  should  provide  the  necessary  boundary 
data.  Appropriate  boundary  values  are  interpolated  from  the  mesh  or  meshes  in  which  the 
boundary  is  embedded.  Thus,  it  is  necessary  to  find  the  interpolation  points  or  cells  on  the 
embedded  mesh  from  which  the  boundary  values  are  interpolated.  Once  these  points  or 
cells  are  found,  their  pointers  are  added  to  a  list  of  such  points  to  be  used  by  the  flow 
solver  to  update  the  variables  at  the  hole  or  boundary  points.  There  are  several  different 
procedures  to  locate  the  points  in  the  embedded  mesh  which  would  be  best  for  interpolation 
boundary  point.  In  general,  the  closet  points  in  the  embedded  mesh  would  provide  the  best 
interpolation  data  for  the  boundary  points,  but  for  some  special  cases,  it  is  necessary  to  pay 
special  attention  in  selecting  the  interpolation  points.  For  viscous  flow  problems,  in  the 
region  very  near  the  surface,  the  closest  points  would  not  necessarily  provide  the  best 
interpation  data,  because  the  flow  properties  change  very  fast  in  the  surface  normal 
direction.  A  small  distance  in  this  direction  would  bring  a  big  change  in  the  flow 
properties.  Therefore,  we  have  to  find  those  points  with  the  same  height  as  the  hole  or 
boundary  points,  and  use  them  as  interpolation  points.  Because  these  points  may  not  be  a 
real  grid  point  in  the  physical  space,  it  is  necessary  to  find  a  method  to  use  the  results  of 
these  real  grid  points  to  get  the  interpolation  data. 

2.2.4  Interpolation  Methods 

Because  the  subdomain  boundary  is  in  the  interior  of  other  computational  domain 
where  the  differential  equations  are  valid,  it  is  appropriate  to  use  the  solution  of  these 
equations  to  update  the  necessary  boundary  data.  There  are  several  approaches  to  obtain 
these  boundary  data,  and  all  of  them  involve  some  form  of  data  interpolation  from  one 
mesh  to  another  mesh. 

We  can  use  tri-linear  interpolation  to  update  the  values  of  boundary  points  for  a 
three-dimensional  problem.  The  tri-hnear  interpolation  which  has  the  form: 
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A  =  tti  +  a24  +  ttjTi  +  +  ttsCri  +  +a8^'nC 

where  0<f,ri,^  <  1  are  the  coordinates  of  the  point  to  be  interpolated  and  the  a^-a^ 
depend  upon  the  function  values  at  the  points  forming  the  interpation  cell. 

In  general,  the  tri-linear  interpolation  can  give  good  interpolated  results.  But  in  the 
region  very  near  the  body  surface  with  very  fine  grids  in  the  normal  direction,  this  method 
may  give  wrong  results  because  of  geometric  and  physical  reasons.  Because  different 
interpolation  methods  need  different  kinds  of  data  structure,  it  is  necessary  to  transfer 
different  kinds  of  data  structure  into  a  standard  data  structure  for  easy  data  management. 

2.2.5  Data  Structure 

A  Chimera  method  requires  the  management  of  a  large  amount  of  grids,  solution 
data  and  interpolation  relationships,  which  consist  of  a  list  of  every  boundary  point  in  each 
grid  which  receives  interpolated  information  from  other  grids,  a  list  of  every  grid  element  in 
each  grid  which  provides  interpolated  information  to  boundary  points,  a  list  of  all  points  in 
the  composite  grid  that  are  hole  and  boundary  (i.e.,  "blanked")  points,  and  the  pointer  from 
boundary  points  to  the  grid  elements  which  provide  the  boundary  point  with  interpolated 
solution  data.  Additional  information,  such  as  grid-storage  location,  number  of  points  in 
each  coordinate  direction,  location  of  interpolation  stencils,  location  of  interpolation 
coefficients,  location  of  interpolated  boundary  data,  and  location  of  hole  or  excluded 
points,  should  also  be  included.  Therefore,  it  is  necessary  to  keep  track  of  the  storage 
locations  of  the  coordinates  of  each  grid,  solution  data  on  each  grid,  all  of  the  interpolation 
information,  and  develop  a  good  management  method  to  transfer  all  these  data  to  the  user  in 
an  automatic  manner. 


21 


Because  most  of  these  data  are  stored  in  lists,  connections  among  the  data  are  made 
through  the  use  of  linked  lists  and  pointer.  It  is  necessary  to  store  the  position  of  the  first 
element  (a  pointer)  in  each  grid  and  the  number  of  points  in  each  coordinate  direction  on 
each  grid  to  locate  the  coordinates  of  any  points. 

2.3  Grid  Generation  and  Composite  Grid 

Although  the  use  of  overset-grid  approach  reduces  the  complexity  of  grid 
generations,  the  task  of  generating  a  satisfactory  body  conforming  grid  about  a  three- 
dimensional  complex  geometric  body  is  still  not  easy.  In  general,  the  grids  must  not  be  too 
distorted.  They  should  have  smooth  variation,  and  should  be  clustered  to  flow  field 
regions,  typically  near  boundary  surfaces.  Moreover,  the  grids  should  be  generated  in  an 
automatic  manner  that  requires  a  minimum  of  user  input 

Because  of  the  application  of  boundary  conditions  can  be  simplified  in  finite 
difference  calculation  if  grid  lines  coincide  with  boundary  lines,  body  conforming 
curvilinear  grids  are  often  used  in  finite  different  flow  field  simulation.  One  approach  for 
generating  body  conforming  grid  is  to  solve  a  set  of  partial  differential  equations  with  a 
series  of  monotone  variation  level  lines  4(x,y,z),  Ti(x,y,z)  and  ^(x,y,z)  as  a  solution  of 
these  partial  differential  equations.  The  values  of  these  variable  ^,r|,and^  are  user- 
specified  on  the  boundary  surface,  and  constraints  expressed  as  differential  equations  are 
used  to  develop  the  grid  away  from  the  boundaries.  The  most  popular  approaches  require 
the  solution  of  a  set  of  elliptic  equations  which  satisfies  the  maximum  principle.  For 
external  aerodynamics  applications,  the  placement  of  the  outer  boundary  curve  will  not 
have  to  be  precisely  specified,  especially  for  Chimera  overset  grid  approach.  Solving  a  set 
of  hyperbolic  or  parabolic  governing  equations  is,  however,  more  powerful  and  efficient. 
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The  governing  grid  generation  equations  are  determined  by  the  differential 
constraints  of  two  mesh  line  orthogonality  relations  between 
^  andri  and  between  t]  and 


and  a  user  specified  mesh  cell  size  constraint 


^(x,y,z) 


=  AV 


where  the  f  is  the  position  vector,  (x,  y,  z).  The  grid  is  obtained  by  first  defining  a  surface 
grid,  then  the  governing  equations  have  to  be  solved  by  marching  outward  in  the  normal 
direction  from  specified  distribution  of  points  on  the  inner  body  boundary,  under  the 
constraint  of  the  two  orthogonal  relations  and  the  user  specified  spacings  or  volumes. 
Details  of  the  hyperbolic  grid  generation  procedure  can  be  found  in  Ref.  55. 

In  the  overset  mesh  technique  a  major  grid  stretches  over  the  entire  field.  The  major 
grid  might  be  a  simple  rectangular  grid  or  a  grid  generated  about  a  dominant  boundary  or 
body  surface.  Minor  grids  are  generated  about  remaining  portions  of  the  body- 
configuration  or  any  other  special  feature  such  as  an  intense  vorticity  region.  The  minor 
grids  are  used  to  resolve  features  of  the  geometry  or  flow  that  are  not  adequately  resolved 
by  the  major  grid.  They  are  generated  somewhat  independently  of  the  major  grid  and  are 
overset  on  top  of  the  major  grid  without  requiring  mesh  boundaries  to  joint  in  any  special 
way.  Overall,  the  effect  of  one  grid  is  imposed  upon  the  other  by  means  of  boundary  data 
that  is  interpolated  back  and  forth  as  the  iterative  solution  process  of  solving  the  difference 
equations  on  each  grid  proceeds.  The  minor  grid  outer  boundary  data  comes  from  the 
major  grid  solution,  not  from  a  far  field  boundary  condition  routine. 
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As  the  separate  grids  are  to  be  treated  as  independent  entities,  boundary  conditions 
must  be  supplied  to  each.  The  boundary  conditions  of  the  differential  equations  which 
model  the  flow  provide  data  only  at  the  boundaries  of  the  computational  domain.  Thus, 
other  data  must  be  obtained  for  the  subdomain  boundaries  which  are  not  coincident  with 
those  of  the  computational  domain.  Because  the  subdomain  boundaries  typically  lie  in  the 
interior  of  the  computational  domain  where  the  differential  equations  are  valid,  it  seems 
appropriate  that  the  solution  of  these  equations,  should  provide  the  necessary  boundary 
data.  Appropriate  boundary  values  are  interpolated  from  the  mesh  or  meshes  in  which  the 
boundary  is  embedded. 

Finding  a  grid  element  which  can  provide  interpolated  information  to  a  hole 
boundary  point  is  a  necessary  but  not  sufficient  condition  for  an  interplation  to  be 
considered  satisfactory.  For  example,  a  grid  which  is  to  supply  interpolated  information  to 
a  boundary  point  may  itself  be  boundary  point  which  is  updated  through  interpolation. 
Most  of  the  information  used  to  update  the  original  boundary  point  will  come  from 
interpolation  rather  than  directly  from  the  flow  solution.  This  close  coupling  of  boundary 
points  can  result  in  poor  global  convergence.  Therefore,  in  the  grid  generation  and 
composite  grids,  a  sufficient  overlap  region  (at  least  two  points)  must  be  provided  and  in 
Chimera  connecting  code,  a  simple  strategy  is  added  to  not  allow  boundary  points  to  be 
updated  by  grid  elements  which  are  themselves  comprised  of  boundary  points. 

2.4  Solution  Procedure 


The  Chimera  approach  has  the  advantage  that  a  flow  simulation  code  developed  for 
a  single  general  curvilinear  grid  can  be  readily  adapted  for  composite  overset  grids.  It  is 
easy  to  write  an  outer  loop  control  program  that  calls  in  grids  and  interface  routines  and 
provides  logic  to  set  communication  relations  between  the  grids.  At  each  time  step,  a  grid 
and  its  data  are  fetched  from  an  isolated  memory  device  into  the  working  memory. 
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Boundary  interface  arrays  that  store  grid  interconnect  data,  .  are  also  brought  in.  The 
Qgf.  arrays  hold  overset-grid  boundary  values  for  the  current  grid  which  are  supplied  from 

the  other  grids.  The  solution  on  the  current  grid  is  then  updated  or  advanced  in  time  using 
the  flow  algorithm,  the  Qgc  interface  boundary  values,  and  the  usual  boundary  condition 

routines  (tangency)  that  are  applicable.  Overset  boundary  data  that  the  current  grid  sends  to 
other  grids  are  then  loaded  into  Qgc  and  all  the  arrays  are  sent  back  to  the  isolated  large 

memory.  The  next  grid  is  then  fetched  and  so  on. 

Because  of  the  above  reasons,  the  intention  is  to  build  a  simple  solver  which  can  be 
used  for  quick  embedment  within  the  Chimera  framework,  check  outputs,  give  good 
guessing,  diagnostics,  etc.,  as  well  as  for  use  with  an  optional  fortified  approach.  These 
numerical  algorithms  for  inviscid  flow  simulation  are  efficient,  and  simple  and  robust  so  as 
to  minimize  the  amount  of  engineering  man  hours  to  obtain  an  initial  solution.  Therefore,  a 
formulation  for  the  steady  inviscid  flow,  in  terms  of  primitive  variables  of  the  original 
Cartesian  velocity  variables  is  described  for  body  fitted  general  curvilinear  coordinates. 
The  details  of  the  formulation  are  given  in  the  later  chapter.  A  limited  number  of 
calculations  have  been  undertaken  to  validate  the  formulation. 
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3.0  FLOW  SOLVER 


3.1  Subsonic  Steady  Inviscid  Nonconservative  Formulation 


Continuity: 


Vorticity  Definitions: 


Ux  +  Vy  +  Wx  =  -■d 


t?  =  {upx  +Vpy+WPx)IP 


=Wy-  Vj, 


U2  —  Uz  —  Wj 


U^3  “■  '^x 


Written  as  a  system  of  first  order  equations 


1  0  0 

0  0  0 

0  0-1 
0  10 


0  1  0 

0  0  1 

0  0  0 

-10  0 


0  0  1 

0-10 
1  0  0 

0  0  0 


Aqx  +  Bq-y  -I-  Cpz  —  g 
/I  0  0  \ 


Multiply  by  D\A*,  with  =  I  0  —1  0  1  then 

Vo  0  -1/ 


1  0  0 
0-10 
0  0-1 


/O  1  0 
+  10  0 
Vo  0  0 


/O  0  1 
+  000 
Vl  0  0 
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or 


M^x  +  B\qy  +  C\qz  —  /i 


Each  matrix.  A\,  Bi.  and  C\  has  real  eigenvalues.  Moreover,  the  matrix  aA\  + 
bB\  +  cC\  has  distinct  real  eigenvalues 


X{aAi  +  bBi  +  cCi )  =  —a,  ±\/ 


and  a  complete  set  of  eigenvectors  must  exist  (Mathematica  could  not  find  them, 
so  they  will  need  to  be  worked  out  )  so  the  formulation  is  applicable  to  general 
coordinates  as  wiU  be  noted  later. 

/-I  0  0\ 

If  Eq.(l)  is  multiplied  by  D2B^,  with  D2  =  I  0  1  0  1  then 

V  0  0  -1/ 


(00l\/ux\  /0  0  0\/«„\  /-I  0  0\  /  Uz\  /-u;2\ 

0  0  oj|n*j  +  |o  0  ®  I 

1  0  0  }  \wx  I  \o  1  Q  I  \wy  J  \0  0  I J  \wz  J  \  J 

(5) 
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For  symmetry,  add  Eqs.(3).  (4)  and  (5) 


/I  1  l\fuA  /-I  1  0\f  Uy\  /-I  0  l\  fUz\ 

1-1  0  +  1  1  1  +  0  -1  1  U, 

\1  0  -Ij  \wj  \0  1  -l)  \WyJ  \l  1  ij  \wj 

((J3  —  CJ2  —  t?  \ 

-  u;3  -  t9  I 

u)2  —  (jJi  —  d  J 

(6) 

or 

Aqr  +  Bqy  +  Cq^  =  f  (7) 

Each  coefficient  matrix  has  real  distinct  eigenvalues  and  the  linear  combination 
aA  +  bB  +  cC  has  the  distinct  real  eigenvalues 

XiaA  +  bB  +  cC)  = -a-b-c,  ±y/Z{a^  +b^  +  c2)  (8) 

and  again  a  complete  set  of  eigenvectors  must  exist. 

3.2  General  Coordinates 

Let  ^  ij  =  T]{x,y,z),  and  C  =  C{x,y,z)  then  Eq.(7)  becomes 

(^x-4  +  ^yB  +  ^eC)q^  +  iVx-A  +  VyB  +  T)xC)qjj  +  (Cz-d  +  (^yB  +  (izC)qQ  =  f  (9) 

All  combinations  of  eigenvalues  follow  from  Eq.(8). 

Rewrite  (9)  as 

Aq^  +  Bqr,  +  Cqc;  =  f  (10) 
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^  _{Up,^Vpr,^Wp^) 

P 

U  =  ^xU  +  (yV-\-^ztV,  V  =  T]:,U-\-r)yV  +  T}^W.  =  Cx^  + 

3.3  Numerical  Algorithm 
Euler  Implicit  {h  =  At  oi  relaxation  factor): 
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(/  +  hA^S^  +  hB^Sr,  +  hC^6^)A^  =  -h[iA^6^  +  8^6^,  +  -  r]  (11) 


where  ^ 

Approximately  factored  in  ADI-like  fashion 

{I  +  hA^6^){I  +  hB^6r,)iI  +  hC^S^)A^  =  -h[{A^S^  +  B^6r,  +  Cr6^)r  -  H  (12) 
Algorithm  form: 

(/  +  hA^6^)A^  =  +  B^6r,  +  C"6c)9"  -  ri 

(/  +  /iB%)A^*  =  Ag* 

{I  +  hC^6(;)A^  =  (13) 

=  ^  +  A^ 

With  smoothing,  above  equations  become: 

(J  +  hA^S^  -  h€i{AV)^)A^  =  -h[(i”A(  +  B^Ar,  +  -  /”] 

-/iee[(AV)2  +  (AV)2  +  (AV)2]g  (14) 

(/  +  hB^6^  -  /i€i(AV)„)Ar*  =  Af 
(I +  hC^8^  -  h€i{AV)^)A^  =  A^* 
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The  left  hand  side  of  first  3  equations  of  (14)  can  be  written  as: 


[/  +  hA^S^  -  /ie,(AV)^]Ag-  =  [{-h^  -  he,).  {I  +  2e,h).  {h^  -  he,)]A^ 


\l  +  hATbr,  -  /ie,(AV)^]  A?-  =  -  he,).  {I  +  2e,h).  (/i|  -  he,)]Ar 


[I  +  hA^S^  -  he,{AV)t^]A^  =  [i-h^  -  he,).{I  +  2eih).{h^  -  he,)]A^ 


3.4  Bernoulli  Equation 


If  the  flow  is  steady,  inviscid,  adiabatic,  the  Bernoulli  equation  can  be  derived 
from  the  Croco  relation  and  perfect  gas  relations  (Anderson  et  al.,  1984): 


I  I  «,f2  \  1  -y— 1 


_  1  4.  ^  (  A/f2  _  ^  ^ 

Qnch  ^  V 


’{9-5oo)IR 


For  irrotational  flow,  the  entropy  correction  term,  e  «oo)/^  can  be  dropped.  To 
avoid  computational  cost  this  is  expanded  using  the  binomial  expansion: 

P  /,  ,  sn  ^  .  n(n-l)  2  ,  n(n-l)(n-2)  , 

=  (1  ±  a)”  =  1  ±  an  +  a^  ±  ^ - ^a^  +  . . . 

Poo  2  2x3 


where 


n  =  1/(7  -  1) 
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Written  as 


Q  = 


7-1  /  ,  ,2  +  V"  + 

—  {M^ - Ti - 


_  (n  —  1) ,  (n  —  2) ,  (n  —  3) 

(1  +  a)”  =  l+an(l+Q^  ^--(1+q^  ,  (1  +  Q  ■-— •  •  •))) 


or 


n(n-l)  _  _^n(n-l)(n-2)  n{n-l)(n-2)(n-3) 

(1+a)"  =  l  +  a(n  +  a( - - - +a( - - - ha - — - ))) 


or 


(1  +a)"  =  1  +a(ci  +a{c2  +  q(c3  +004))) 


where 


Cl  = 


7-1 


Therefore 


Poo 


^  1  \ 

C2  =  iri - r  -  1) 

2  7  —  1 

I  ^  o\ 

<^3=  3(:^-2) 

CH  =  f(-i-r-3) 

4  7  -  1 


=  1  +  a(ci  +  a(c2  +  a(c3  +  004))) 


3.5  Boundary  Conditions 


Boundary  conditions  are  taken  from  a  vorticity/tangency  inviscid  relation. 
From  the  vorticity  definitions: 

Cy^C  -  CzVC  =9l  =(z}i-  i^yW^  +  TJyWr,  -  “  VzVr,) 
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(zH  -  Cx«’<:  =  92=^2  -  {izH  +  VzU-^  -  -  VxUJrj) 

Ci(-’c  ~  Cy^c  =  93  =  ^2  -  +  nxVr)  -  iyU(_  -  riyU-T)) 

On  the  wall,  these  equations  can  be  written  as 

Cy«’  -  C,zV  =  -AC^l  +  Cy  f^’  -  =  Pi 

Gu  -  C,xW  =  -ACg2  +  CzU*  -  C,xW*  =  92 

C,xV  -  C,yU  =  -ACP3  +  CiV*  -  CyW*  =  P3 

where  u*.u*  and  lo*  are  values  of  u.v,  and  w  at  AC  above  the  wall. 

The  condition  on  the  surface  C  =  o  is  given  by 

CiU  “1“  Cy^  ^ 

Solving  these  four  equations  for  the  three  unknowns  u.  ik  and  w  on  the  surface  by 
generalized  inverse  (or  using  the  vorticity  relations  to  remove  the  other  components 
from  the  tangency  relation)  gives  the  vorticity/tangency  relations: 

(Cl  +  Cy  +  Cf  =  CzP2  —  CyP3 

(Cj  +  Cj+<l)f  =  C.93-C.9i 

id  +  Cy  +  C.)'»  =  Cv9i  -  Cih 

or  in  vector  form 

(Cx  +  Cy+Cz)9  =  -VCxJ 
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4.0  VERIFICATION  AND  APPLICATION 

4.1  Potential  Solution  for  Flow  around  a  Sphere 

The  ideal  flow  over  a  sphere  is  given  by  the  superposition  of  a  uniform  stream 
and  a  doublet.  The  stream  function  of  these  contributions  can  be  written  as: 

0  =  [(^)  -Q 

in  which  tq  is  radius  of  the  sphere,  and  it  is  easy  to  verify  that  for  arbitrary  9. 
r  =  ro  is  a  spherical  stream  surface  where  4’  =  0-  The  velocity  components  for  this 
flow  are: 


The  maximum  occurs  at  the  equator,  where  q  =  —3Uj2.  The  surface  pressures  are 
found  by  Bernoulli’s  equation  on  the  surface  of  the  sphere,  in  terms  of  the  pressure 


coefficient,  the  result  is 
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We  can  transfer  the  solutions  (V^.  V#)  in  spherical  coordinate  system  to  (u. 
of  Cartesian  coordinate  system  with  following  relations: 


r  =  \/x^  -h  -h 


cosO  = 


sind  =  v/l  —  cos^O 


costp  = 


sinil>  = 


rsinO 


then: 


u  =  VrCosd  +  Vffsintp 


V  =  VrSinBcosi)  +  V^cosOcostp 


w  =  Vrsin\dsintl)  +  VecosQsin'il’ 


With  these  relations,  for  any  given  grid  point  position  (x,y,z),  it  is  easy  to  have  the 
potential  solution  u,v,  and  w. 


4.2  Solutions  for  Flow  around  a  Sphere  in  Single  Grid 


Because  the  analytical  solution  for  a  potential  flow  around  a  sphere  is  available 
for  comparison,  and  the  sphere  grid  is  relatively  easy  to  generate,  the  flow  pass 
a  sphere  was  selected  as  the  model  to  verify  the  simple  inviscid  flow  solver  given 


above.  Theoretically,  the  flow  solver  given  above  should  give  the  same  result  as 
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potential  flow  solution  as  long  as  the  vorticity  terms  disappear  from  the  source  terms  in 
equation  (1).  The  Mach  number  contours  on  sphere  surface  and  a  cross  section  in 
computational  space  are  given  in  Figure  10  through  13.  All  of  these  results  show  good 
agreement  with  the  analytical  results,  therefore,  it  was  concluded  that  the  codes  were 
functioning  properly. 
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4.3  Flow  Solutions  with  Multioverset  Grids 


Since  the  simple  inviscid  flow  solver  and  connecting  codes  are  working,  it  is 
necessary  to  test  these  codes  by  calculating  relative  complex  geometric  configurations.  We 
have  a  hyperbolic  surface  grid  generation  program  and  the  3-D  hyperbolic  grid  generation 
program.  For  a  complex  multibody  configuration,  it  is  not  too  difficult  to  use  these  grid 
generators  to  generate  3-D  grid  on  each  body,  and  superimpose  them  together,  and  use  the 
Chimera  method  to  solve  the  whole  problem.  Even  just  using  the  simple  inviscid  flow 
solver,  it  still  can  get  some  good  approximate  results  which  would  be  very  useful  for 
further  complex  and  detailed  calculations,  such  as  solving  boundary  layer  or  Navier  Stokes 
equations.  It  is  not  difficult  to  introduce  different  flow  solvers  for  different  domains  or 
meshes,  which  could  treat  very  complex  problems  both  in  geometry  and  flow  properties. 
This  is  an  important  advantage  of  Chimera  method. 

Figures  14  and  15  show  the  Mach  number  contours  of  the  solution  of  a  flow  over  a 
sphere  covered  by  four  pieces  of  grid,  at  two  different  Mach  numbers.  One  is  at  M„  =  0.2 
and  the  other  is  at  its  critical  Mach  number  M.  =  0.55.  Both  figures  show  good 
agreement  with  potential  solutions. 
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Figure  14.  Mach  Number  Contour  on  Sphere  Surface  for  four  grids  at  M 
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Figure  15.  Mach  Number  Contour  on  Sphere  Surface  Interpolated  from  Four  Grids 
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4.4  Flow  Solutions  with  Two-Body  Configurations 


In  general,  for  complex  multibody  configurations  problems,  there  will  be  some 
holes  in  the  mesh  when  the  solid  bodies  are  embedded  in  the  mesh.  The  simples  3-D  case 
is  two  spheres  next  to  each  other,  and  each  grid  is  embedded  in  the  other  grid,  as  shown  in 
a  two-dimensional  cross  section  in  Figure  16.  There  is  a  hole  in  each  mesh  where  the  other 
sphere  bodies  are  embedded.  Figure  17  shows  the  velocity  vectors  in  a  symmetric  section 
of  the  flow  field  for  the  two  spheres.  Two  spheres  and  ellipsoids  placed  next  to  each  other 
are  also  calculated  as  a  test  case  at  =0.20,  which  is  in  Figure  18  through  20.  Both 
results  are  reasonable,  therefore,  it  is  possible  to  say  that  the  code  can  handle  multibody 
configuration  problems. 
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Figure  1C:  One  of  two  grids  with  a  hole  in  it  for  two  sphere 
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Figure  18;  Pressure  (;()uto\ir  in  a  cross  section  of  two  spheres 


MACH  NUMHER 
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Figure  19:  Pressure  contour  in  a  cross  section  of  two  spheres 


CONTOU 
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5.0  CONCLUDING  REMARKS 


A  complete  computer  code  has  been  described  that  implement  3-D  grid-embedding 
techniques  as  a  part  of  a  flexible  solution  concept  that  we  have  called  a  Chimera  method. 
The  code  utilize  procedures  for  combining  grids,  locating  embedded  boundaries  and 
interpolation  points,  and  manipulating  complex  data  structures.  The  validity  of  the  method 
was  successfully  demonstrated  on  several  geometries  for  inviscid  flow. 

A  simple  and  robust  algorithm  for  steady  inviscid  flow,  in  terms  of  primitive 
variables,  is  described  for  body  fitted  coordinates.  This  algorithm  has  been  adapted  for  the 
Chimera  framework  for  the  calculation  of  complex  geometric  problems.  The  flow  about  a 
sphere  and  an  ellipsoid,  for  both  single  and  overset  grids,  showed  good  agreement  with  the 
exact  solution  for  subcritical  and  critical  cases.  Comparison  with  an  implicit  Euler  solver 
for  the  flow  about  an  ellipsoid  showed  the  computational  advantages  of  the  present  method. 
Therefore,  the  flow  solver  presented  here  could  possibly  be  quite  useful  for  preliminary 
engineering  analysis. 
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